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Introduction 

An understanding of chemically reactive flows is necessary in a wide range of disciplines 
including astrophysics, plasma physics and the chemical industry. Recently, applications in 
aerospace engineering have lead to a number of studies aimed at deriving phenomenological equa- 
tions, the Navier-Stokes equations coupled to reactions, from the Boltzmann equation for increas- 
ingly complex systems including internal degrees of freedom and three-body interactions[l, 2]. 
However, these investigations do not exhaust the range of interesting applications. A number of 
important applications arise in the physio-chemistry of cavitating bubbles[3]. Aside from obvious 
examples such as the study of flames and explosions, a recent area of interest is sonochemistry in 
which ultrasound is used to induce conditions of extreme temperature and pressure inside bubbles 
with the effect of dramatically increasing the rates of chemical reactions[4]. Closely related is 
the phenomena of somnoluminescence - in which a fluid irradiated with ultrasound emits light - 
is believed to be caused by pressure waves acting on small bubbles of gas in the fluid(see, e.g. 
[5]). The bubbles are subjected to such rapid compression that Shockwaves may develop and the 
concentrated energy drives many chemical reactions particularly when the shocks reach the center 
of the bubbles giving rise to high temperatures and densities. In fact, it has been suggested[6] that 
some (endothermic) reactions may play an important role in limiting the temperatures reached in 
the center of the bubble. It is therefore of interest for these applications, as well as some of the 
others mentioned above, to understand the phenomenological equations governing a reacting gas 
under extreme conditions and far from equilibrium. The Boltzmann equation cannot be considered 
an adequate basis for such a study due to the fact that it is only applicable at asymptotically low 
densities. In fact, the only simple fluid that is amenable to analytic investigation at finite densities 
is one composed of hard spheres. The purpose of this paper is therefore to review the kinetic the- 
ory of reacting hard-sphere systems and to use this as a basis for a hydrodynamic description of a 
reacting fluid far from equilibrium. In particular, the kinetic theory will be used to derive the exact 
balance equations describing the local concentration, density, velocity and temperature fields from 
which the extension of the Navier-Stokes equations to include reactions is developed based on the 
Chapman-Enskog procedure applied to the Enskog approximation to the kinetic theory. A primary 
result will be to show that the usual phenomenological description consisting of the Navier-Stokes 
equations coupled to a reaction-diffusion-advection equation is only applicable if the chemical re- 
actions take place on a time scale which is comparable to the dissipative time scale Xk"^, where A 
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is a transport coefficient and k a typical wavevector. If the reactions are slower, then all hydrody- 
namic relaxation takes place before the chemistry gets started and chemistry and hydrodynamics 
are effectively decoupled. Faster reactions, with time-scale comparable to ck, where c is the sound 
velocity, leads to additional couplings of the reactions to the hydrodynamic fields. Even faster 
reactions lead to the chemistry taking place so fast that hydrodynamics is irrelevant. 

The hard-sphere interaction model has proven remarkably useful as models of single- and 
multi-component simple fluids since, in many respects, the phenomenology of the hard-sphere 
systems and atoms interacting via more realistic pair potentials is qualitatively identical. For ex- 
ample, hard-sphere systems exhibit the full range of transport coefficients found in all simple 
fiuids[7], possess a freezing transition[8] and the structure of hard- sphere fluids in equilibrium is 
not much different from that of any other fluid[9]. On the theoretical side, the equation of state 
of hard-sphere fluids is easily modeled[9], kinetic theory is simplified by the fact that only binary 
collisions are important and it is possible to formulate an extension of the Boltzmann equation - 
the so-called Revised Enskog Theory or RET [10]- which not only describes the transport prop- 
erties of multi-component hard-sphere fluids at finite densities but which also describes transport 
in the solid state[ll]. More recently, inelastic hard spheres have been used as a model for driven 
granular fluids with similar success. The hard-sphere interaction is therefore an ideal model for 
understanding the extreme conditions occuring in sonochemical experiments. 

The kinetic theory of chemically reacting hard spheres has in fact been discussed in the 
literature[12-15]. The principal aim of these studies was to investigate contributions to the re- 
action rates coming from dense fluid effects (e.g., ring kinetic theory leading to mode-coupling 
models) at equilibrium. In these studies, the atoms carry labels indicating their species (some- 
times called their color) and all intrinsic properties like the atomic radius and mass is specific to 
the species. When the atoms collide, there is a probability that a reaction takes place in which 
the species labels, and hence atomic properties, change. The probability typically depends on the 
rest-frame energy of the colliding atoms: if the rest frame energy is greater than some specified 
activation energy, the reaction can take place with a probability that is, in general, a function of 
the relative energy. Energy may be gained or lost (exo- or endothermic reactions) but the sizes of 
the atoms are generally invariant since, were they to also vary, a collision could result in one of 
the atoms overlapping with third atom. (Technically, there is no reason that atoms could not get 
smaller and most results would apply to such a model). Besides being restricted to a chemistry 
consisting of color labels (and so excluding, e.g., the exchange of mass upon collision), a common 
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assumption in earlier work is that the chemical reactions are slow compared to other transport pro- 
cesses. Since the rate of chemical reactions is generally determined by the ratio of the temperature 
to the activation energies and the difference in the concentrations, this implies that the results are 

only applicable near equilibrium or for low temperatures. One of the primary goals of the present 
work is to indicate how the phenomenological description (Navier-Stokes equations coupled to 
a set of advective-reaction-diffusion equations) must be modified to account for large deviations 
from equilibrium. 

The organization of this paper is as follows. The second Section develops the formal statistical 
description of the system. The discussion of possible collision rules is followed by the develop- 
ment of the Liouville equation and the exact balance laws describing the evolution of the local 
mass, energy, momentum and partial densities. The Enskog approximation is also introduced. The 
third Section discusses the Chapman-Enskog solution of the Enskog equation and, particularly, the 
difference between the assumptions of fast and slow chemical reactions. It is shown that, if the 
chemical reactions are sufficiently slow, the fluid may be described by the Navier-Stokes equations 
for the total mass, energy and momentum densities and a reaction-diffusion-convection equation 
for the concentrations with the only coupling between the two being the convective term occur- 
ring in the latter (i.e., the "usual" description). However, under less restrictive assumptions, the 
reactions are shown to depend in a much more complicated way on the hydrodynamic fields. The 
paper concludes with a discussion of the physical meaning of the different assumptions. 

I. STATISTICAL MECHANICS OF REACTING HARD-SPHERES 

Consider a system of hard spheres of various species confined to a volume V with positions 
{~qi}f^i and momenta {'^i}f^i- Each atom will also be described by a set of discrete labels, k 
for the i-th atom, which fix its "chemical" properties (mass, hard-sphere diameter and reaction 
parameters). When two atoms collide, both the mechanical variables and the species labels of the 
atoms change. The dynamical variables are altered according to some deterministic collision rule 
so that for collisions between the i-th and j— th atoms, 

^ — ^ijj^i (1) 

~^ -^j ~ 
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where the coUision operator bij. describes a coUision involving the reaction li + Ij — > + I'j 
(i.e., the i-th atom changes from species k to l[, etc.) and clearly, one expects that h^J, = h^.^,. 
Because the species change instantaneously and randomly upon collision, the species labels must 
be viewed as discrete random variables. Attention here will be restricted to the model in which the 

probability of making a particular transition is given by some function of the relative phases of the 

I'V- 

two colliding atoms K{^^, (xij) where the notation indicates that this probability depends on the 
relative velocity, Ij'ij — Ij'i — Ij' j, and position ~qij (e.g., through the combination (Ij'ij ■ ~qijf). 



A. Collision Rules 



The coupling between chemistry and hydrodynamics can be captured simply by considering 
atoms that carry a label (e.g., color) that can change during collisions. Allowing for the non- 
conservation of energy gives a relatively broad model that includes endo- and exothermic reac- 
tions. However, in the interests of generality, the problem of modeling reactions that not only 
violate energy conservation, but that also allow for the exchange, or even loss, of mass will be 
considered. 

The modelling of the collision rules in the case that mass is either exchanged or lost upon 
collision is somewhat problematic. To understand why, consider the usual arguments leading to 
specular collision rules in the case that mass is invariant. Defining the total and relative momenta 
respectively as 

P = ^1 + ^2 (2) 
= P^l - 'P2, 

the conservation of total momentum means that 

P' ^ P (3) 

P ^ P + 1 llll 

where 7 ^^^^ (^1,2^2) = — 7 ilil {^2, xi) is to be determined. Second, the energy balance equation 
can be written as 

E{x[ , 4) + SE^l^^ (xi ,X2)^ E{x[ , 4) (4) 

I' I' 

where SEi^i'^ {xi, X2) is the energy lost during the collision. Substitution of Eq.(3) gives 



where the reduced mass is u;,/, = "^h'^h -dimensions, this gives one constraint on the D 
independent components of 7 ^^^^ so, for example, it fixes the magnitude of 7 ^^^^ if its direction is 
known: In one dimension, the problem is therefore solved. In higher dimensions, the conservation 
of angular momentum gives the needed additional constraint. This reads 

P' X 2l±l^ + 1^;, X tl2 = P X ll±ll + 1^, X tl2 (6) 

or, using the conservation of total momentum, 

^11 X ti2 = (7) 

^i' I' 

thus fixing the direction of 7 ^^^^ as being along the line joining the centers of the atoms. In this 
case, eq.(5) gives 



TilZ = "^f^hh { -^12 • "g 12 - y (^12 • ~qi2f - -j^SEllllj qi2. (8) 

I' I' 

Taking SEi^^^ — gives the usual result for elastic hard spheres whereas setting the loss to be a 

fixed fraction of the contribution to the rest-frame kinetic energy due to the velocity along the line 

II II II II ^ ^ 2 

joining the atoms, = \^i^^'hh ( ^ 12 ■ <1 12) , is the model used for inelastic hard spheres 

v v I v v 

(i.e., granular fluids) and gives a coefficient of restitution a^^j^^ = W 1 — A^^^^ 

When mass can be exchanged upon collision, it is useful to introduce the total mass M/^/j = 
fnu + fni„, the center of mass Q — "^'^ ^A^""'^ ^ ^ and the center of mass velocity V — P /Mi i 
Notice that even with the conservation of total mass, the center of mass is not generally invariant 
if mass is exchanged and the positions are kept fixed. So, the conservation of angular momentum 
gives 

Mi'^i'^V' xQ' + ixii^iiv' X ^12 = Mi^i^V xQ+ iii^i^v x ^12 (9) 
or, using the conservation of total momentum, 

Mi^i^V X {q' - q) + (/x/i^^^' - 111,1,-^) X = 0. (10) 

In general, this equation cannot be satisfied since it implies 

'vx fg'-g)VVi2 = o (11) 



which is not generally true. The conclusion is that any collision rule which does not leave invariant 
the center of mass will necessarily result in a violation of the conservation of angular momentum 



at the microscopic level. It is not possible to compensate by allowing the positions of the atoms 
to shift during the collision since this could lead to overlapping configurations involving a third 
atom. In fact, one would expect that the inclusion of internal degrees of freedom, in particular of 
rotation of the spherical atoms, would allow for a galilean-invariant collision law and this will be 
explored at a later date. For present purposes, given that the collision rule cannot be uniquely fixed 
by appeaUng to general principles, the only recourse is to try to construct reasonable models. One 
possibility is to conserve the angular momentum in the center of mass (CM) rest frame since then 
V — and angular momentum can indeed be conserved. Another is to work in analogy to the 

case of invariant masses and to require that either all momentum transfer be along ~qi2 (so 'p'' — 

^ I' I' ^ 

p + 2fj,i-^i,^^i^i^qi2) However, since the former is not Galilean invariant when mass is transferred 
(since in a frame moving at velocity li the relative momentum is 'p' boosted = + — ^h) ~^)^ 
it is not clear how to uniquely apply it. In fact, if one tries to enforce this constraint in the CM 
frame, it gives the same result as fixing the angular momentum in the CM frame. A second 
possibility is to demand that all velocity change be along q 12 (so v' — v + ^ililqi2). For 
illustrative purposes, both options will be considered. 

For the sake of generality, it is also interesting to consider the consequences when mass is not 
only exchanged, but is lost. It is clear that the previous considerations concerning the specification 
of the collision rule under mass exchange apply to this case as well so that, again, a model must 
be introduced in order to specify the relation between the lost mass and the total momentum and 
energy. Without considering specific applications, it is not clear that any unique conclusions can 
be drawn, so by way of illustration, I will assume that the mass is carried away in such a way 
that the total momentum in the CM frame is conserved. This means in general that the law of 
conservation of momentum becomes 

P' = P - 5mllv (12) 

I' I' 

where ^m^J^^ = mi^ + mi^ — mi'^ — mi'^. The model adopted here is that the mass is carried away 

by n^J,^ particles with masses m° so that (^m^J^^ = X]i=i^ with rest-frame velocities 

/A 

satisfying Xlj=i^ — 0- the case in which the angular momentum in the rest CM frame is 

held constant, it is natural to also require that these particles carry no net angular momentum. 

I' I' 

Finally, some model must be specified for the energy lost (or gained), SEi^i^, which might in- 
clude contributions due to kinetic energy that is carried away by the lost mass and energy lost (or 
gained) through other mechanisms (excitation of internal degrees of freedom, radiation, exother- 
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mic and endothermic chemical reactions...)- Suppressing the species indices for a moment, the 
energy differential can be written as the sum of two contributions, SE — SE^n + SEq, where 
the first is the energy carried away with the lost mass and the second is due to any other in- 



elastic processes. In the CM rest frame, 5Em = XliLi 1"^* C^'i) = so in the lab frame, 
5Em = Yl^=i ("^i* ^ ^) ^ ^^"^ ^ \5mV'^. Further, I assume, as is commonly done, that 
the remaining energy loss (or gain) is frame independent (which means in particular that it can 
only be a function of the relative velocity of the colliding atoms). The energy balance equation 
therefore reads 

^ -Pi + + + -^mftv^ = + (13) 



where it is understood that ^E^^^^ — bE^^^^{xx^x-^ is a galilean-invariant function of the phases. 
This expression depends on the model for the lost mass (if any) but is independent of any other 
assumptions concerning the collision rule. 

From Eq.(lO), conservation of angular momentum in the rest frame then gives 



12 — ,, 



(^12 + A;i;^g) (14) 



for some scalar {xi, X2). Substituting into Eq.(13) gives 



>S: J • Q)" - (1 - V- - - ^Smllv^ (15) 

and one then finds that 

1111 f mv — mu mv — mu \ — >• 

= 2^i,iXt^,, + \ - \ ]p- (16) 

Demanding that the change in the relative velocity be along the line joining the atoms gives a very 
similar result 

;' V 



\2 = ^12 + Xitllm (17) 



with 



= -^12 ■ gi2 



;' /' 



2 Sni 

\ K^i2 ■ mf - (1 - ^^i.iJ^i'a) vl, - -—SEfl^ - -^V^ (18) 



and 



= _ 1 + 2«,,,A;;i;f:, + - '-Vi ) I P (19) 
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Note that these two models coincide in the special case that the reduced mass is invariant which 

obtains in one of two circumstances: the atomic masses are invariant or if the atoms just exchange 

I' I' 

masses so that m;^ = mi^ and vice-versa. In both cases, mass is necessarily conserved, ^m^^^^ = 0, 
and all other conclusions are model-independent consequences of galilean-invariance. 

B. The evolution of phase functions 

The dynamics of any hard-sphere model consists of free streaming interrupted by binary colli- 
sions. In non-reacting fluids, the collisions lead to an instantaneous change of the velocities of the 
colliding atoms. The generalization to the reacting fluid only requires that the chemical species 
labels, and hence the masses and any other species-specific properties, to be viewed as dynamical 
variables as well and, so, as part of an enlarged phase space. 

Two atoms, say atoms i and j, collide at time Tij when their centers are separated by their 
relative hard- sphere diameter al^l. 

I ~Qi {^ij ) - ~qj {jij ) I = (Tjj ) = all. (20) 

where, e.g., for additive models, the relative hard sphere diameter is simply the sum of the atoms' 
radii u;.;^ = + \(yij- The atoms do not have to all have the same size (e.g., an acceptable 
possibility is that different species have different sizes but chemical reactions always transform 
atoms from a species of a given size to other species of the same size). An exceptional possibility 
in which size could change is one in which atoms only get smaller upon collision: this might be 
useful to model certain granular materials that fragment upon collision (e.g., the ice composing 
the rings of Saturn) and could be handled within the present formalism as long as the position of 
the center of mass of each atom is invariant. From Eq.(20), one has that the time of collision is 




where the sign has been chosen according to give the physical solution. If the right hand side is 
imaginary, then no collision takes place for the given initial conditions. This aspect of the dynam- 
ics is independent of what actually happens after the collision and is the reason that the structure of 
the pseudo-Liouville equation is independent of the collision rule. The pseudo-Liouville equation 
describing the time evolution of an arbitrary phase function, A (F; = A (xi, /i, ...xjv, In'-, t) then 
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follows immediately by analogy with the non-reacting fluid and is 

= I^A + C,A (22) 

where the binary collision operators are 

^+ (u) = - Vij • 'vijS {qij - ai.i.) e {-Qij ■ 'Vij) M^^ - 1 . (23) 



'J 



As discussed in Appendix A, this can be derived directly for a system of two atoms by writing the 

exact solution to the two-body problem and differentiating; the generalization to A^-atoms follows 

I'l'- 

immediately due to the fact that only binary collisions occur. Here, M.i'.i. (xij) is a random matrix 
which, in any realization, takes on the value 1 for some single combination of l[, Ij and is zero 
otherwise and which is distributed according to 

(-Mm:(^..)) =K^i^^^)■ (24) 

(The notation used here indicates stochastic quantities by means of Calligraphic type and uses 
carets to denote operators and averages over the stochastic process are denoted as {■■■),.^acr^ ^ 

V-V- 

non-reacting system, it becomes M-H {xij) = 6i'_i-6i'j^ . The only other formal difference from the 
non-reacting case is that the momentum transfer operator, bi./., has the effect of altering both the 
mechanical variables and the species labels. So, just as this operator instantaneously changes the 
position in phase space of the ith atom from before a collision at time t to Xi{t+) — x'^{t^) 

it also instantaneously alters the species labels from k to k = the difference being 
that x[{t-) is a deterministic function of Xi{t-) and Xj{t-) whereas the evolution of /-(t-) is 
stochastic. For phase functions which have no explicit time dependence, the Liouville equation 
can be formally solved to get 

^(r, t) = exp (c+t^ A{r) (25) 

which has the meaning that the system evolves from the initial phase T. 

The most important difference from the non-reacting system appears in the evaluation of statis- 
tical averages. In the presence of reactions there are two statistical processes that must be consid- 
ered: the distribution of initial conditions and the stochastic process that alters species labels at the 
collisions. For a given distribution of initial conditions (F) = p\^j^ (xi,X2, ...xn) (giving 
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the probability that the first atom begins with species h and phase xi, and so), one has 

{A; t) = jdv (r) ^ jdv (r) (exp (Z^t) ^(r))_^^ (26) 

and the notation should be understood as implying a sum over the initial species labels 

/ drp(°)(r)/exp(£+t)A(r)) (27) 

J \ V / / react 

= ^ J dxi...dxN pf^lj^ {xi, ...Xn) (exp (^C+t^^ ^^^^A{r) 

h.-.l]^ 

where A(r) can be taken outside of the average over reactions since it depends only on the ini- 
tial conditions. Now, since each collision involves an independent stochastic process, it follows 
that /cxp [jC^tj \ = exp t) which is evaluated using Eq.(24). Thus, the time 

\ V / / react V \ / react ' 

averages become 

(A; t)^ J d^i-d^N pfl,i^ {xu -Xn) exp (l+t) A{T) (28) 

h...lN 

with the deterministic operator 



' dxi 



and the (reaction-averaged) collision operators are 



This shows that, from the point of view of evaluating the statistical averages, it suffices to work 
with the deterministic dynamics defined by L+ which no longer treats the species labels as discrete 
stochastic variables. Instead, the phase functions are at all times averaged over the reactions and 

so do not explicitly represent dynamical quantities as might be realized in a computer simulation. 
In fact, they correspond to the average result of an ensemble of simulations, all beginning with 

I'-V- 

identical initial conditions, but differing in the realization of the reaction process M.^ {xij). 



C. The evolution of the distribution function 

The adjoint of the Liouville operator L+ is defined as 

J d^B (r) L+A (r) = j dr (l^b (r)) a (r) (3i) 
11 



from which one finds (see appendix B) 



i 

with the adjoint collision operator 



(32) 



i<j 



with 



. o,,b 



e {-l^ij ■ ~qij) S {qij - ai.i.) 'Vij ■ 'q^j 

(33) 



(34) 



Here, the operator ( fe^y 1 is the inverse of 6^^,^ both in terms of the change of the mechanical 
variables as well as the species labels so that for an arbitrary function (V^l^^ B {xi, If, xj, Ij) — 



loll 



B y\hah j ct; \^ahj ■ illustrate the structure of this operator, consider the case 

of inelastic hard spheres used to model granular fluids. Specializing to a single species, one has 



^ij = ^^ij = - (1 + «) '^ij ■ Qij^ 



where a is a constant, from which it follows that 



givmg 



d b Xi. b Xj 



1 + a 



a 



1 + a 



a 



^ ij ■ Qij 



1 

a 



(35) 



(36) 



(37) 



so that 

(ij) B {xiJi;Xj,lj) 



-b-' - 1 
a 



e {-l^ij ■ ~qij) 6 {qij - ai-i.) 'Vij ■ ~qijB (xj, k; Xj, Ij) . 

(38) 

which is the usual result[16]. 

An important generalization of this result concerns the case that the inverse transformation 
V^l^Xi is not unique. This can happen even in the single species, inelastic case if the coefficient of 
restitution depends on the velocities. For example, if a = a {^ij ■ qij) then the inverse collision 
rule is determined by solving 



^ ij ' ^ij o; ( f ■ qij) V ij • qij 



(39) 
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which may or may not have a unique solution. In the latter case, must be written in terms 
of a sum over the various branches and must include step functions which restrict the domain of 
integration in Eq.(31) to the appropriate domain for each branch. In practical calculations, it is 
usually most convenient to recast integrals over (ij) into integrals involving T+ (ij) so as to 
avoid this complication. 

Given the adjoint operator L^, Eq.(28) can be written as 



{A; t) = j dV (exp {lH) (F)) A{T) = j dV p (F; t) A (F) 



(40) 



where the second equality defines the time-dependent distribution function. It's time dependence 
is given by the pseudo-Liouville equation 

where in the standard notation[7] 

T_ (ij) = -T^ (ij) . (42) 

The Bom-Bogoliubov-Green-Kirkwood-Yvon (BBGKY) hierarchy follows immediately from 
the Liouville equation. Defining the reduced distribution functions as 



fh...lrn {Xl...Xm) 



— > i dxm+i-.dxN p{T) . (43) 

(A^ - m)\ ^ J 

lm+l---lN 



and integrating the pseudo-Liouville equation over Xm+i---XN and summing over the correspond- 
ing species labels gives the m-th equation of the hierarchy 



^ + E^^-^+ E T.{zj)]fi,..,^{x,...xJ (44) 

, 1=1 ^ * l<i<j<m / 

m „ 

= - E E / ^^rn+iT- {im + 1) fl^...lm+i {Xi...X^+i) . 

The first equation of the hierarchy is the starting point for the Enskog kinetic theory as described 
below. 



n. EXACT BALANCE EQUATIONS 



Now consider the phenomenology of the reacting fluid which is expressed in terms of the 
macroscopic hydrodynamic fields. The results presented here are derived using only the general 
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form of the collision rule, eq.(3), and the microscopic energy balance equation, eq.(13) so that 
the only assumptions made with respect to the collision model are those concerning the energy 
transported by any lost mass. 

The local fields of interest are the number fractions 



ni ( r", t) = 6 (T" - 'qi) 5^,; = J (Pui fi 



fi{T,ir,;t), (45) 

and the mass, momentum and energy densities, defined respectively as 

p{'r,t) = ^m,n,(r',t) (46) 

^n{'r',t)kBT{'r,t) = ^mz^l/^^^ ( - Vi) ; = X] j dl/iV^fi{V,lfi;t) 

where the excess velocity is Vi{t) — l^i{t) — li {~qi, t)and the total number density is 

I 

It is also convenient to introduce the number fractions, or concentrations, xi{~r,t) — 
ni t) /n (7^, t). The balance equations for these quantities follow directly from their defi- 
nitions and the first equation of the BBGKY hierarchy. The details of the derivation are given in 
appendix C and only the results summarized here. 

A. Number, mass and concentration 

Integrating over the positions and velocities gives the balance equation for the local partial 
number density 

+ V • {-Ttm) + V-7i^ S}""^ (48) 

with the source 

^f"-* ( r", ^) = J d^id^2 {~qi2 ■ ^12) S {qi2 - (Ji^i^) e {-qi2 ■ ^12) (49) 

ablil2 

xfhh {xi, X2; t)5{'r' - ~qi) K^i^ {X12) {5ai + Sbi - Su^ - Su^) 
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and the number current j i — j f + j Y with 

7f t)^ J dir, fi (r, ir,,t) (50) 

and 

7r (7", t) = -^Yl J dxidx2 ~qi2 i~qi2 ■ ^12) 5 (^12 - (^hh) © (-^12 ■ ^12) (51) 

ablil2 

xfhh {xi,X2; t) K^^i^ {X12) {Sal - Sbi - Sii^ + 611^) 

X dx 6 (7^ — x~qi — {1 — x) ~q2) ■ 
Jo 

The source term represents the gain or loss of atoms of type / due to chemical reactions. The kinetic 
part of the number current is familiar from the study of multi-component, non-reacting systems[7] 
where it takes the form j f — DijVrij + LiVT + o (V^) and, e.g., gives rise to Pick's law 
when substituted into eq.(48) . Here, it is seen that this diffusive current is enhanced by a second 
contribution, eq.(51), that arises solely from the reactions (i.e., it vanishes if K^-^i^ — Sai^Sbir^). 
This is due to the transport of type-i atoms due to the reaction process. The conservation of total 
number density immediately follows by summing over the species label 

V • (^n) + V • J]7f = (52) 

I 

where the sum over the species of the coUisional contributions to the number current vanishes. 
Similarly, multiplying by mi and then summing gives the balance equation for local mass density 

^^p+V ■{lip) + V -Q ^S^f\ (53) 

where, the mass flux is 

Q ( r", ^) = XI / ^^i^^2 ~qi2 {~qi2 ■ ^12) S {qi2 - CThh) (-$12 • ^12) (54) 

abhh 

xfhh {xi:X2; t) {X12) {ma -mb- mi^ + m^J 

X / dx S — x~qi — {1 — x) ~q2) , 
Jo 

which vanishes if no mass is transported during collisions, and the mass source term is 

Sj-''\V,t) = ^ X] y dxidx2 {~qi2-^i2)S{qi2- aij^i^)Q{-qi2-l^i2) (55) 

ablil2 

xfhh {^1, X2\ t)5{'r - -qi) K^X (^12) ^TT^hh^ 
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which is only nonzero if the coUisions do not conserve mass. Finally, using the definition of the 
concentrations, xi — ni/n, the reaction equation is found to be 

d 



dt 



xi-\-li ■ V xi + n ^ 



■ j i-xiV - Y, 3 



K 



loin) 



where the term on the right is now identified as the reaction rate. 



(56) 



B. Momentum and velocity fields 



The balance equation for the local momentum, written in terms of the local velocity, is 

—pit + V • {pltli) + V-\^P+Qltj = S + It si"' 



(57) 



with the pressure tensor P — P ^ + P^ + P ^ where the kinetic contribution is 

P'^(y,i) = ^m, f V^iV^i, (58) 

I 

and the collisional contribution is 

P ^ , ^) = XI / ^^i^^2 'qvi {'qi2 ■ ^12) 5 (gi2 - © {-qi2 • ^12) (59) 

xfhh {xi, X2; t) kIiIi {xu) / da; 5 ( y - x~qi - {1 - x) ^2) 



, 7' /' 

where 7 ^^^^ is the change of momentum in the rest frame 

~lhil = + (^i'l - - m'^ - mi^) Vi2- (60) 
Finally, the contribution from the instantaneous exchange of mass is 

-P^ = "2 XI / ^^1^^2 "?^12 Cg"l2 • ^12) (?12 - O-ilia) Q (-$12 • ^12) (61) 

lll2l[l'2 

xfhh {X1X2) k\^iI {X12) (Vi2-li^ {mv^ - mi^ - mv^ - mi^) 

X dx 5 ("r^ — x~qi — {1 — x) ~q2) 
Jo 

The source terms arise due to momentum being carried away by the lost mass and the new term is 

given by 

^ip) _^ 1 ^ f _^ ^ _^ 

S ( ' ^) ^ 2 / ^^^^^2 ( ? 12 • vi2)S {qi2 - (Ti^h) © (-gi2 • V 12) (62) 

X {Vi2 - ^) fhh {xi,X2] t) kIiIi {X12) Sm\l\lS (t - ~qi) . 
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By using the balance equation for total mass density the equation of motion for the velocity field 
is found to be 

^'u + li ■ Vli + (v - ^ +Q ■ V^) = p-^^ . (63) 



C. Energy density and temperature 

The balance equation for total energy density is 

where the new source term, arising if energy is not conserved by the collisions, is 

= ^ / dxxdxi Cq 12 ■ ^ 12) S {qi2 - (Ti^i2)Q{-qi2- ^12) (65) 

hhl[l'2 



X 



;/ 7/ 1 7/ 7/ / — > , \ ' 



fhh {xi,X2; t) kIiI^ (X12) S{y - 'qi) , 



which is recognized as the generalization of the source term studied in the context of granular 
fluids. The heat flux is written as a sum of several contributions 

where the kinetic part has the usual form 

I 

as does the first part of the coUisional contribution 

'q^ (y , ^) = XI / '^^i^^s 'q\2 Cq 12 • ^12) S {qu - ai^i^) Q {-qu • ^12) (67) 

hl2l[l'2 

xfhh {xi,X2; t) K^l^l (X12) [y - uy i^i^ dxS{r - x q i - {1 - x) q 2) , 

which is a measure of energy displacement during the collision (i.e., one atom experiences a net 
gain of energy, the other a net loss and this represents an instantaneous movement of energy from 
the location of the second atom to the location of the first). Qualitatively new contributions arise 
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from the instantaneous transfer of mass. 



^1 mi>mi^-mi>mi^ 

X j dxidx2 ~qi2 {'qi2 ■ ^12) 5 {qi2 - ai^i^) 6 (-$12 • ^12) 

Xfhh {XUX2] t) Kilil {X12) Hi^i^v 

X dx 5 ("r* — x~qi — {1 — x) ~q2) , 
io 



and from the loss of energy, 



1 mi[ - mi'^ 

2 ^ my + my 

lll2l[l'2 1 2 



X J dxidx2 ~qi2 Cq 12 ■ ^12) S {qi2 - ai^i^) © {-qi2 ■ ^12) 



^fhh {xi,X2] t) Kilil (X12) 6Eilll 

X / dx 6 {1^ — x~qi — {1 — x)~q2) ■ 
Jo 

Alternatively, noting the relation between the total energy and the kinetic temperature 

the evolution of the kinetic temperature is found to be given by 

: + V • "g" 



dt 



I 



Dnk} 



Dnki 



(68) 



(69) 



(70) 



(71) 



D. Enskog Approximation 

The expressions for the balance equations are exact. As a consequence, they depend on both 

the exact one-body and two-body distribution functions which are, in principle, determined by the 
BBGKY hierarchy. For example the equation for the one-body distribution is explicitly 



d ^ d 
at q 



a,b,l2 



V 2 



(72) 



X 



d{Xi,X2) 

X© (-^12 • V12) S {qi2 - CTjijJ ^12 • $12 
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However, since the latter cannot be solved exactly, except in the special case of equilibrium, it is 
necessary to introduce an approximation. The most common approximation is to assume that the 
velocities of two colliding atoms are uncorrected prior to the collision (they are of course corre- 
lated after the collision since the collision itself generates correlations). That this approximation 
is sufficient to decouple the BBGKY hierarchy is seen from the fact that the right hand side of eq. 
(72) since the step function, 6 (—^12 • V 12), is non-zero only for atoms approaching one another 
and the delta function restricts the domain to the instant of contact. Thus, the assumption that 
atoms are uncorrected just prior to a collision, Boltzmann's "assumption of molecular chaos", is 
precisely the statement that 

e (-^12 • 'qi2) 5 (?12 - (Tl^h) fhh CQI: ^1, 2, ^2) (73) 

~ e {-~Ui2 ■ 12) 6 {qi2 - ai^i^) g {~qi, 2; t) fi^ {xi; t) fi^ (2:2; t) 

which, when substituted into Eq.(72) gives the Enskog approximation to the one-body distribution 
function. The factor of g {~qi, V2; t), the spatial pair distribution function, allows for spatial cor- 
relations which always exist. In the Revised Enskog Theory, it is approximated by the equilibrium 
functional of the density evaluated for the local density field of the fluid[10]. The same approxi- 
mation can be above to give the corresponding Enskog-approximation to the balance equations. A 
final consequence follows from the second equation of the BBGKY hierarchy which has the form 

= -njd3 {T _{13) + T _{23)) f 1,1,1, {xi,X2,X3;t) . 

Since atoms cannot interpenetrate the two-body distribution must have the form fi-^.^ (xi, ^2; t) = 
6 (gi2 - CT/^/J Vhh i^iy ^2] t) for some function yi^i, {xi, X2; t) which is continuous at qi2 = ai^i,. 
Then, one expects that the singular terms in Eq.(74), arising from the gradient acting on the step 
function and from the definition of T_(12), must cancel, gives the constraint 

'Vl2 ■ {qi2 - CTl^h) fhh {Xl: t) = -T_{12)Ji,i, (xi, X2\ t) (75) 

and some rearrangement, together with the approximation of Eq.(74) gives 

5 {qi2 - (Jhh) fhh {xi,X2; t) ~ 5 (gi2 - ai^i,) g ("g i, "g 2; t) fi, {xi, t) fi^ (^2; t) (76) 

- (^12 • $12)"^ T_{12)g ("g 1, 'g'2; t) fi, {xi;t) fi, {X2; t) 
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which expresses the two-body distribution function at contact in terms of a completely uncorre- 
lated piece, the first term on the right, and a correction that takes into account velocity correlations 
generated by the collision, the second term on the right. This can be used to evaluate two-body 

correlations at the Enskog level of approximation [17-19]. 

m. CHAPMAN-ENSKOG SOLUTION 

In the previous Section, the exact balance equations were developed and the Enskog approx- 
imation introduced. Next, this framework is used to derive the explicit equations governing the 
evolution of the hydrodynamic fields by means of the Chapman-Enskog approximation. As noted 
in the Introduction, previous studies of the kinetic theory for reacting systems have often made 

the assumption that the chemical reactions are slow relative to the hydrodynamic time-scales. The 
primary goal here is draw out, and make more precise, the meaning of this condition by outlin- 
ing the Chapman-Enskog procedure under different assumptions about the speed of the chemical 
reactions. 

Before beginning, note that the phrase "hydrodynamic fields" usually refers to those local fields 
which are conserved in the long-wavelength limit (which is to say that their sum over the entire 
system is conserved). For a non-reacting fluid of hard-spheres, this means the local partial number 
densities, and the momentum and energy densities. For a reacting fluid, the partial number den- 
sities are not conserved and for models of endo-/exothermic reactions, even the energy is not be 
conserved. Following the practice developed in the study of granular fluids (which are non-reactive 
but do not conserve energy) it seems natural to expand the definition of "hydrodynamic" fields to 
include those fields which would be conserved in the limit of vanishing reaction probabilities. A 
partial justification for this is that all of these fields are necessary to develop a meaningful descrip- 
tion of the non-reacting fluid, so one expects that they must also be included in any description of 
the reacting fluid (i.e., a minimal-coupling argument based on continuity of the description with 
respect to the control parameters). 

The Chapman-Enskog procedure attempts to construct a so-called normal solution of the En- 
skog equation which is to say, a solution which is a local functional of the (exact) hydrodynamic 
fields and for which all of the space and time dependence occurs implicitly through those fields[7]. 
This implies that the space and time derivatives of the distribution function can be written in terms 
of the corresponding derivatives of the fields and the functional derivative of the distribution with 
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respect to the fields. In other words, one has 

fa CQi, 'v'l,t) = fa [vilXi (Vl, t) , 71 {~qi, t) , ll ^Qi, t) , T (^1, t)] 



(77) 



so all of the dependence on g i and t occur through the hydrodynamic fields so that the time 
derivative can be expressed as 



d 

Q^fa{~qi,^l,t)^Yl 



dxi dfa ^ dndfa ^ &u_ _ dfa ^ dT dfa 



(78) 



dt dxi dt On dt du dt dT 
Then, the kinetic equation determines the functional dependence of the distribution on the fields 
and their derivatives, while the fields are in turn fixed self-consistently by the balance equations 
(in the Enskog approximation). 

A further approximation which is made in practical calculations is to assume that spatial gradi- 
ents are small so that the equations can be solved perturbatively via a gradient expansion. To order 
the terms, one introduces a uniformity parameter e and replace V with e V and order terms in e. 
Since the space and time derivatives are related by the balance equations, one also introduces an 
expansion of the time derivative = dt = d^^^ + ed^^^ + ... as well as of the distribution itself 



fa {~qi,l^i,t) = fl ['vi\xi,n,li,T\ + efl [v i\xi,n,li ,T\ + 



(79) 



where the notation indicates that the distribution is a functional of the hydrodynamic fields. These 
expansions are substituted into both the Enskog equation and the balance equations and an order- 
by-order solution is sought. Writing the Enskog equation as 



d d \ 

^ + V^i • fa {Xi]t) = Jab,cd [/c, fd] 

^^1/ bed 



(80) 



so that 

Jab,cd \fci fd\ 



I 



d 



(b''J,xi,blix2 



bit) Kt - Sac6> 



bd 



d{xi,X2) 

xQ {-~uu ■ T12) S {qu - (Ted) ^12 ■ qi2fc (Ti, ^1; t) fd (T2, ^2; t) 



(81) 



it is also necessary to expand the non-locality of the collision operator which comes from the term 
5 {qi2 — cTcd) — ^ (512) + ^cFcd^' (?i2) + ■■■where the derivative of the delta-function, which will 
give rise to spatial gradients of the distribution, are scaled with an appropriate factor of e. In order 
to control the speed of the chemistry relative to the hydrodynamics, the non-diagonal part of the 
reaction probabilities is separated out as 



Kf, ^ 5a.5bd + - ^achd) 



(82) 
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giving 

Jr -f -f 1 JC JC j{invariat) r r i* 1 i en j{veactive) r r /■ i /0'3\ 
ab,cd [Jc, Jd\ — OacObdJafi U«' Jb\ + ^ -'a6,cd Uc, JdJ (.OJj 

where the non-reactive, or invariant, part is the usual collision operator for non-reactive (but pos- 
sible energy non-conserving) multi-component fluids 

Ja,b [fa: fb\ = d q 2d V 2 

XO (-^12 • 12) S {qi2 - (Tab) 'Vl2 • Qllfa (Yl, ^li ^) fb {~q 2, '^2, t) 

and, as indicated, the reactive part of the collision operator will be arbitrarily treated as being of 
order a in the gradient expansion. Thus, the full expansion of the collision operator will take the 
form 

T ^ £ £ ^ r r / Y{invariant)0 r <■ /• 1 , j(invariant)l r j- j- ^ , \ /acs 

Jab,cd[fc, fd] ^ OacObdyJa,b [/a, /ft] + eJ^b [fa,fb] + -j (85) 

_L^a f Ti'^eactive)0 rr f 1 , ^ j{reactive)l rr f 1 , \ 
+^ ['^ab,cd i/o Jd\ + ^Jab,cd i/o Jd\ + -) ■ 

A. Zeroth-order 

The zeroth order equation for the distribution is then 

qO^(O) j{elasUc)0 r rQ /.Q] , r j(reacUve)0 r rQ /.Q] /q/;\ 

"tJa - 2^ 'Ja,b [Ja > /b J + ^^aO 2^ J ab,cd [Jc ' Jd\ (-OO) 

b bed 

which must be supplemented by the corresponding equations for the fields expanded to zeroth 
order 

d^xi = 5,ori-'S'("^(''^"^**"^)° (87) 
d^n = 
d^Tt = 

QOrj-^ ^{invariant)0 _j_ ^ ^^{reactive)0 

These balance equations, together with the assumption of normality, eq. (78), serve to define the 
meaning of the term in eq.(86). Note that the fluxes do not enter, being of first order in the 

gradients, and that the sources are separated into a non-reactive and reactive part using eq.(82). For 
the concentration, mass and velocity fields, there are in general no non-reactive contributions to the 
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d(b:ix,,b:ix2 

d{Xi,X2) 



(84) 



sources whereas for the temperature, there is the possibiUty of such a contribution in which case, 
one recovers the inelastic hard-sphere system used to model granular fluids. Furthermore, use has 

^zf{p){reactive)0 

been made of the fact that /i ' must be a function of — 1? which implies that S = 

(since there are no other zeroth-order vectors) so that no source can appear in the velocity equation, 
at this order. These zeroth-order equations illustrate a complication that occurs for fast reacting 
systems (e.g. a — 0) compared to non-reacting multi-component systems: namely that the sources 
in the balance equations at order n require knowledge of the n-th order distribution. For non- 
reacting elastic systems, the n-th order balance equations generally require only the n — 1 order 
distribution so that there is no coupling between the two. Non-reacting inelastic systems, i.e. 
granular fluids, share this complication as can be seen from the appearance of the source term 

^{non)0 in Eqs.(87). 

For a > 0, only the temperature can have a zeroth-order time dependence and so can contribute 
to the left side of eq.(86). If this temperature source is zero then the left hand side of eq.(86) is zero 
and the fa^^ will simply be proportional to a Maxwellian. The solution of eqs.(86) and (87) for the 
case a > and the non-reactive source in the temperature equation being non-zero corresponds to 
the so-called homogeneous cooling state in granular fluids and has been discussed in detail in the 
literature for single-component[20] and multi-component systems[21]. 



B. First-order 



At first order, one has 

d^fi'^ + {dl + ^1 • V) /° (88) 

El j{invariant)0 r nQ /■!] , j{invariant)0 r nl /■O] _i_ j{invariant)l r eO xOl \ 

b 

I r \ ^ / j{reactive)0 r pQ nil . j{reactive)0 r nl nQl . j{reactive)l r nO fO~\\ 

+OaO y-'ab,cd [Jc ' /dj + '^ab,cd [jc ' /dj + '^ab,cd [Jc ^ Jd\ ) 

bed 

,r j(reacUve)0 r rQ /.Q] 



bed 
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and for the fields 



djn + V • 



un) 



d}T + li ■VT + 



Dnkf 



dill + li ■ Vli + p ^Vp 



(0) 



r —1 ri{n){reactive)l 
. r —1 Q{n){reactive)0 



_^^(p)(o) 

5aOP S 
^ (invariant) 1 

I ^ ^^(reactive)! _j_ ^ ^^(reactive)0 



(89) 



where we have used the fact that at zeroth order there are no velocity-independent vectors and 
only the unit tensor available so that we must have P = p^'^^ 1 , W — w^'^^ 1 and all vector 
fluxes must vanish. 

In general, the first order distribution must take the form 



(r,V;t)^nxi<j>i (V 



hi{V)+Ai{V)V ■ \/n + BiV ■ VT 

+ C", : (yV - ^Tv • + (90) 

where I have written the zeroth order distribution in the form //°^ = nxi(j)i (V). Here the coeffi- 
cients Ai,Bi, ...are scalar functions of the velocity (and in general depend also on space and time 
through a dependence on the local hydrodynamic variables as does the zeroth-order distribution, 
although for the sake of conciseness this dependence has been suppressed). The function h{V) 
represents the first-order correction to {V) due to the energy-dependent chemical reactions: for 
example, if the only allowed interaction were A + A ^ A + B and this only took place if the CM 
kinetic energy were greater than some threshold, Eab, then starting with a system of all ^-type 
atoms, one would expect to build up a preponderance of fast B atoms and a corresponding deficit 
of fast A atoms. It vanishes in the case that a > 2 and energy is conserved by the non-reactive 
dynamics. The consequences of different orderings of the reaction terms will be considered sepa- 
rately. 



1. Ultra-slow reactions: a > 2 



In this case, there are no reactive terms in the first order equations. The solution is therefore 
the same as for the equilibrium (or HCS) multi-component system. The second order balance 
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equations will also have no reactive terms. Summing up to second order, the Navier-Stokes order 
balance equations are then 



d_ 
di 



xi + li ■ Vxi + n ^ 



d 



V • 7"; -XiV ■ ^ j 



K 







(91) 



— n + V • + V • J] J f = 



+ 'u ■ Vli + f V • - ^ V • Q 

ot V 







dt 



y : + V • "g + V • (^^ • ty) + • Vu^ 

(92) 



where the fluxes are the sum of zeroth- and first-order contributions, P — P + P and 
the source ^ consists of (non-reactive) contributions summed through second order. This means 
that for the granular case, ^(*"^««^««"^*) ^ Q, the Navier-Stokes order balance equations require 
knowledge of the second-order (or Burnett order) distribution function. There is, at this order, no 
coupling between the hydrodynamic equations and the reaction equations. Inserted into Eq.(91), 
the result would be the Navier-Stokes equations for elastic hard-spheres, or their generalization 
for inelastic hard spheres. If this expansion is continued, the a — th order balance equation for the 
concentrations would be 



■K{a-l) 
I 



— 1 ci{n){reactive)0 



(93) 



Clearly, the reaction equation remains unknown in this case since one would need to consis- 
tently include the higher order hydrodynamic contributions that would come from the number 
current which, in turn, would bring in couplings to higher-order gradients of the hydrodynamic 
fields. Without knowledge of these higher order terms (and they are not known for even the 
one-component fluid) the reaction equation can only be consistently studied in the absence of 
hydrodynamic gradients when the reactive terms are treated as of order a > 2. 
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2. Slow reactions: a = 2 



In this case, the first order solution is again the same as in the non-reacting case. However, the 
sources will have second-order contributions so that the Navier-Stokes equations take the form 



—xi + u ■ Vxi + n 
at 



_ ^-1 g{n){reactive)0 ^^^^ 



+ V • + V • J] 7f = 

I 

Q ^ , ^ ^ ^ ^ ^^(ip)(reactive)0 

—It + It ■ Vlt + p-^ (V ■ P - itV ■ Qj ^p-^S 



T + ^ • VT - - V • y 7f + — ^ : + V ■~q + ^ ■ (it -w] + 1-Q ■ Vu^ 
n ^ Dukb I \ / 2 _ 

^{invariant) _j_ ^{reactive)0 (95) 



For the simplest case that the reactions conserve energy and momentum, the reactions are governed 
by exactly the convective-reaction-diffusion equation that one might expect. The reaction rates are 
calculated using the local equilibrium distribution as in elementary treatments [22]. Except for 
the usual modification of the transport properties arising from the use of the Enskog equation, as 
opposed to the Boltzmann equation, there are no new dense-fluid effects. 

3. Moderate reactions: a = 1 

For moderately fast reactions, the situation becomes more interesting. Considering here only 
the case that mass and energy are conserved by all collisions, the first order balance equations - 
the generalization of the Euler equations - are found to be 

at ) 

d -> _ 
—n + V ■ un = 
at 

(^^ + lt -V^lt + p-^Vp^^^ = 

so that the reactions, with reaction rates calculated from the local-equilibrium distribution func- 
tion, enter into the Euler equations. The Navier-Stokes equations will involve the reaction rates 
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calculated up to first order in the distribution. In general, the only nonzero coupling in the reac- 
tion source will take the form _ 5'^") y • li where 5";^"^ is a scalar function of the 
concentrations, density and temperature. The Navier-Stokes equations will therefore take the form 

I + ^ • X, + || V • 7, - X, V • E 7f ) = ^-lg(n)ireacUve)0 ^^^^ 

_|_ ^-1 ^{nKreactive)! _|_ 

+ V • + V • J] 7f = 
(^^ + lt -^^11 + p-^V • = 

The source term for the reactions has three contributions: the zeroth order reaction rate (calculated 
using the local equilibrium distribution function), the first order correction (due to deviations of 
the distribution from local equilibrium) and a new, dense-fluid effect which couples the reactions 
to the divergence of the velocity field with some field-dependent coefficient, 5*;*^"^ . This coupling 
is a dense-fluid effect which does not exist in the Boltzmann approximation and not surprisingly, 
its origin is closely related to that of the bulk- viscosity which is also zero in the Boltzmann theory, 
but not the Enskog theory. (The calculation of these terms will be discussed in detail in a future 
publication but the fact that these are the only possible couplings is due to the fact that no other 
galilean-invariant scalars, linear in the gradients of the fields, can be constructed). 

4. Fast reactions: a = 

In the case of fast reactions, no a priori assumption is made about the speed of the reactions 
compared to the hydrodynamic time scales. The balance equations to first order, i.e. the Euler 
equations, are 

— +U-V ]n + nV -u = 

ot 



^ + ir- V ]T + 



dt 7 Dnk 



py . ^ _|_ y . (^w'u) — ^(*"''"''*""*)0 _|_ ^{reactive)0 _|_ ^{tnvartant)! _|_ ^ 
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so that even the Euler equations show the dense-fluid correction to the reaction rates. The second- 
order, or Navier-Stokes, equations require the evaluation of the source terms to second order, which 
in turn requires knowledge of the distribution function to second order (also called Burnett order). 
One then expects that even for mass and energy conserving interactions, the reaction equation will 
contain couplings to gradients of all of the hydrodynamic fields. However, since the complete 
Burnett-order Chapman-Enskog solution of the Enskog equation is not even known for the case of 
a single-component fluid, there is no practical value in continuing the analysis for this case. 

IV. CONCLUSION 

In this paper, the kinetic theory of reactive hard-core systems has been extended to include 
the possibility of mass transfer and/or loss and energy gain/loss. When mass is not conserved, 
the collision rule becomes dependent on the model used to describe the lost mass. Nevertheless, 
quite general expressions for the dynamics of phase functions and the distribution function of the 
system can be given and used to derive equally general expressions for the exact balance laws for 
mass, momentum, energy and concentration. For example, when mass is conserved but a fixed 
fraction of the rest-frame kinetic energy is lost during collisions, the usual inelastic hard sphere 
kinetic theory, used as a model of granular fluids [16], is recovered. The kinetic theory was used, 
within the Enskog approximation, to discuss the various phenomenological laws, extensions of 
the Navier-Stokes equations, that arise from different orderings of the reaction terms within the 
Chapman-Enskog procedure. It was noted that the intuitive model of the Navier-Stokes equations 
coupled to a reaction-diffusion equation through a convective term only arises when the "speed" 
of the chemical reactions is comparable to some hydrodynamic time scale and even in this case, 
an additional coupling to the divergence of the velocity field can occur in dense fluids. 

How fast do we expect the chemistry to be relative to the hydrodynamics? In the model con- 
sidered here, chemical reactions cannot be faster than the collision time. In fact, a typical reaction 
rate would be something like p {5x) e~^^'"^'^Ucoi where p is the probability of a reaction occuring 
if the colliding atoms have energy greater than the reaction energy barrier, E, 5xis the difference 
between the concentration of the species and its equilibrium concentration and u^oi is the collision 
frequency. On the other hand, the Chapman-Enskog procedure is based on a gradient expansion, 
the small parameter e will generally be a measure of the ratio of the typical microscopic length 
scale, the mean free path Imfp, to a typical length scale for hydrodynamic gradients L. (In Fourier 
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space, where gradients V correspond to wavevectors k , this becomes e ~ kl^fp-) So, setting 

p^-E/ksT^^^ ~ {klmfpT J^col ~ (Imfp/L)" V^ol giveS 

Q, ^ _L^J (99) 

ln(Z^/p/L) 

For systems in which hydrodynamics is appUcable, one has Imfp/L << 1 so that a ranges from 
a minumum of > for fcoT ^ E, to very large values for low temperatures. Far 

ln(i,„/p/-Lj 

from chemical equilibrium, 5a; ~ 1, the lower limit could be arbitrarily close to zero depending on 
the reaction probability so that "moderate" and "fast" reactions are possible at high temperatures, 
hideed, if all of these parameters are fixed, then fast reactions will always occur in the hydodry- 
namic regime limit Imfp/L — > 0. The conclusion is that, unless the concentrations are close to 
their equilibrium values, the reaction probabilities are very small or the temperature is extremely 
low, the concept of slow chemical reactions may be of limited applicability and, so, the correct 
phenomenological description, from the standpoint of kinetic theory, may be more complex than 
the reaction-diffusion-advection model. 

In summary,if chemical reactions are slow compared to the rate of dissipation in the fluid, then 
hydrodynamics and chemistry are not meaningfully coupled. If the reaction rate is comparable to 
the rate of dissipation in the fluid, i.e. AA;^ for some transport coefficient A and wavevector k, then 
the usual reaction-diffusion-advection equation results. For faster reactions, additional couplings 
occur and the chemistry and hydrodynamics become more interdependent. The detailed solution 
of the Enskog, and the resulting phenomenological equations for particular reaction models will 
be the subject of a future publication where the importance for sonochemistry of additional terms, 
such as those occurring in Eq.(97), will be investigated. 
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APPENDIX A: THE HARD-CORE LIOUVILLE OPERATOR 



29 



The goal in this appendix is to provide motivation for the statement in the text that the form of 
the pseudo-Liouville operator is independent of the coUision rule. To start with, restrict attention 
to a system of 2 atoms. Let X {T;t) be the characteristic function for collisions after at time t 

beginning with the phase T at time so that if the two atoms do not collide between during the 
interval [0, t] then X (F; t) = whereas if they do collide, X (F; t) = 1. Then the time evolution 
of the phase function Ar (t) — A (F(i), is given by 

{t) = (1 - X (F; t j) A (Fo(t), + X (F; t) A (F'(t), t) (Al) 

where Fo(t) is just the phase of the system propagated a time t into the future in the absence of in- 
teractions and is explicitly Fo(t) = + ~vit^ 'v'l, li, ~q2 + ^2^, ^2, ^2^ • (Note that attention 
is restricted to the case that velocities are constant during free-streaming: generalization to include 

one-body forces is straightforward.) The phase point F'(t) is the position the system would reach 
in phase space if a collision occurred at some time r (F) G [0, t]. Explicit expressions can also be 
given for its components such as ~q[{t) — ~qi + 'v'iT + {t — r), etc. Direct differentiation 
then gives 

dAr (t) _dA dq{t) d 

+ '^-^^^iM^'(tU)-A{Vo{t\t)). 
Now, from the definitions above 



(1 - X (r; t))Y,v,- ^Sj-A (r„(f), i) + X (r; t)J2^:- ^SjrrA (r'(f), 

_^ A (J 



j=l,2 



givmg 



Now, since X (F; t) has the form of a step function (it is zero if t < r (F) and one otherwise), we 
must have 

dX{r;t) 



dt 
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^5(t-r(r)), (A5) 



and this also gives the correct result (zero) if r (F) is imaginary (indicating that no collision ever 
occurs starting from the given state). Then, using 



6{t-r (r)) {A (r'itlt) - A {ro{t),t)) ^ 6{t-r (r)) {A (r'(r (r)), r (r)) - A (ro(r (r)), r m) 

^5{t-r (r)) I 5] Mf£y;^ -UArir (F)) 



v v 



gives 



dAr {t) 



9i 



1=1,2 



;' /' 



(t) (A7) 



I order to express the right hand side entirely in terms of r{t) rather than the initial condition F 
the temporal delta-function is rewritten using 



S(qi2 (t)-ai,i^) =^ 



^912 (^)|t=^.(r) 



(A8) 



where (F) are the roots of i'^i) ~ ~ 0- They correspond to the time at which the two 
atoms are first in contact, i.e. the physical collision time which is denoted as r (F), and the time 
at which they are last in contact if they are allowed to pass through one another (which is not 
physical). One picks out the correct root by noting that at the physical collision time ~qi2 (r) • 
^12(7") < while the sign is reversed at the unphysical collision time so that 

N^/ ^ /X -> .NX (Ht-T(T)) (Ut-T(T)) ()(f-T(T)) 

I ft^i2 [t) I ^^^(r) I ('t) • 1^ 12 (r) I I qi2 {t) ■ v 12 {t) \ 



or 



5(t-T (F)) = S (gi2 (t) - ai,i,) e (-^12 (t) ■ ^i2(t)) \qi2 (t) ■ ^i2(t)| 
giving finally 



with 



dAr (t) 
Jt 



l<i<2 



'd_ 
di 



Ar{t) 



d 



+ E ^+(12; ^) 



(A9) 



(AlO) 



(All) 



(A12) 



l<i<j<2 
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where, for arbitrary phase function B (T^t), 

f+{12;t)B{r{t),t) = S{qu{t)-ai,i,)e{-ti2{t)-l^i2{t))\qi2{t)-l^i2{t)\ (A13) 



X 



For more than two atoms, one simply extends the sums in eq.(A12) since, in a finite system, only 
binary collisions can occur. 

Starting with an initial condition Ar (t) — A (F), iteration of eq.(Al 1) immediately gives 



= VA (F) (A14) 



with L = L(0) which implies that 

{t) = exp {L?j A (F) (A15) 

and 

d 



—At (t) = L exp (Zt^ A (F) = LAr (t) (A16) 



as claimed in the text. 



APPENDIX B: THE ADJOINT LIOUVILLE OPERATOR 

To derive the adjoint operator, begin with its definition 

J dV B{T)L+A{T) = J dV {L^B{T))A{T) (B1) 

or, more explicitly, 

J dxidx2... B{r)L+A{r) = J dxidx2 {L^B{r))A{r) (B2) 

li,l2--- h,h 

Now, 

L+ = Lf + Y,T+itj) (B3) 

i<j 

and it is obvious that, neglecting surface terms, 

JdFB (F) LfA (F) ^ Jdr (^-L^B (F)) A (F) (B4) 
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so 



(0) 



(B5) 



Next, consider one of the collision operators and restrict attention to a system of two atoms. 
Then 



/ dXidX2 B {xi,li]X2,l2)\T+{12)A{xi,li]X2,l2)] 

= - ^ I dxidx2 B {xi, h; X2, k)^ {-^12 ■~q 12) HQi2 - CTu) If' 12- ~q 12 

h,l2,a,b ' 

X {^hh (^1' ^2, k) A$1^i^xi, a;bf^i^X2, b) - SaiM^A {xi, h; X2, h)^ 
Consider the first term. Relabeling the species in the sum gives 



(B6) 



- 5Z y dXidX2 B (Xi, h] X2, I2) 6 (-^12 ■ '^'12) 5 {qi2 - CThh) ^12 ■ 'qi2 
h,l2,a,b 

xK^X (^1' ^2, k) Aihf^i^x,, a-htli^X2, h) (B7) 
= - ^ J dxidx2 B {xi, a; X2, b)Q{-l^ 12- 'q 12)3 {qi2- CTab) 1^12 ■~q 12 

xK':^^ (xi, a; X2. b) A^^j^^x,, liM'x2, h). 

Assuming that the collision operator is invertible, then introducing new integration variables = 
^ab^^i corresponding Jacobian 



ab 



y2 



d{yi,y2) 



(B8) 
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gives 



- ^ y dxidx2 B (xi, h] X2, h) © (-^12 • V12) S {qi2 - ai^i^) l^i2 ■ 'q 

h,h,a,h 



12 



(B9) 



i2\ \._.i,]q( _(yih^ 



ab 1/2; 



Vl2- qi2 



xS iqi2-crab) (b 



'ab 



12 • g'i2i^^,'^ fe) " yi, a; fe) " 1/2, 6 Z^; y^, ^2) 



- 5Z dyidy2 A{yi,li;y2,l2) 

J^;'^ (yi, y2) (?.^,'^) ~' © (-^12 • ti2) S {q,2 - ai,i,) • ~qi2Kj^b' (yi, h; y2, ^i) 

{xi,h;X2,l2) 

where the operator (V'^^'^^ has the effect of changing the species from li, I2 to a, h. One can then 



write 



Ti{l2)B{x^,h;x2M) 



(BIO) 



= - E W (^1' ^2) i^^^ (a;i, /i; 0:2, /2) - 1 

ab 



© (-^12 • 'g'12) 5 (gi2 - C^iiia) ^12 • V 



12 



XS {xi,li\X2,l2) 



In some cases of interest, the coUision dynamics may not be invertible. For example, suppose 
that coUisions with total rest frame energy less than some threshold, E, are elastic while those 
with energy greater than this are inelastic. Then, a pair of atoms with rest frame energy after 
collision of |/Xi2^'i2 < ^ might have resulted from either (a) a collision between two atoms with 
rest-frame energy below the threshold or (b) a collision between two atoms that had energy above 
the threshold but that lost part of this due to the inelastic process. In this case, it is necessary in 
the definition of the adjoint operator to include an additional sum over the various branches of 
the inverse collision dynamics. Even when it occurs, such a complication may not be of practical 
importance since it is often the case that expressions involving the adjoint operator can be 
rewritten in terms of the original operator T+. 
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APPENDIX C: DERIVATION OF THE BALANCE EQUATIONS 

In this section, the general form of the local balance equations is derive and then specialized to 
the partial density, momentum and energy fields. 

1. General form of the balance equations 

Consider any one-body phase function of the form 

(^) = E ^■'Pk (^^) s {r - ti) hi (ci) 

i 

and its average 

= n y" d'vifi{'r,^i,t)ipi{'vi) 
The balance equation for this follows from the first BBGKY equation 

(^^ + l?i • (xi) = E / (12) (0:1x2) (C3) 

and is 

(r, t) + V • / d^i (r, ^1, t) ^1^/;; (^1) (C4) 

= ^ / dXi (^1) (j^ - / dX2T _ {12) f 11^ {x-^Xi) . 

Introducing the specific velocity, V li/r ,t) = ^1— (V,t), the second term on the left becomes 

j dPvi fi{'r,^i,t)l^iipi{'vi) (C5) 

while it proves more convenient to rewrite the right hand side in terms of the r+ collision operator 

/ dxi Sii^i/Ji^ {~vi) 6{r - ~qi) / dx2T_ (12) fu^ {X1X2) (C6) 

hh 

^ X] / dxidx2fii^ {X1X2) T+ (12) ( v^i) 5 (r - 'qi) 
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so that the balance equation becomes 

{V,t) + V -It {V, t) *z {V,t) + W ■ J (Tvi fi {V, -vi, t) ViiJi (^i) (C7) 
^ ^ J ^^^^^^fhh ixiX2) S{r- i) T+ (12) dii^ipi^ {'v'l) 



d 
~dt 



with 



T+{12)6u,'4'h{^i) (C8) 
-~qi2 ■ (qi2 - (Thh) © (-$12 • ^12) ^ K'^il {X12) h}^^^ - 1 5^,-0^1 (^1) 



7' P 



= -"g 12 • 'vi25 {qi2 - cThh) © (-$12 • ^12) XI {^hh (^12) b\l\l - Si^i'^Si^i'^^ Su.ijji^ (v'l) . 

'1*2 

In general, the right hand side can be separated into a sum of a flux and a source term. Let 

Bi,i,,v^,^.,i {x„ X2) = (<g (0:12) bll - Si,i,Si,,^) 5u,i^i, (^1) (C9) 
and define its even and odd components as 

Fhh,i'A^ (^1' ^2) = \ {Bhi2,i'A;i i^ux-i) - Bi.i^^i'^i'^.i {x2, xi)) (CIO) 
= \ {^hh (^12) - Si^i'^Si^i'^^ {6ii,iJi, (^1) - Sii.i^i^ (^2)) 

Shh,l[l'^;l {Xl, X2) = ^ {Bl,l^,l'^l'^-l (Xi, X2) + Bi^i^^i'^i'^.i (X2, Xi)) 



[k^ (X12) tl^t - 



so that 

d 



^^^i{T,t) + V ■u{T,t)^iiT,t) (Cll) 
+ V ■ J] /" dlTi fi (r, ^1, t) ViiPi (^1) 

- X / dxidx2fu2 {X1X2) (5 ( - g'l) "g'12 • ^"12^ (gi2 - ai^i^) 6 (-$12 • ^12) 
X (^iii2,riii;i (^1,2:2) + Si^i^^p^i^^-i {xi,X2)) . 
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Then, relabel the dummy variables to give 

jdxidx2fi^i^{xiX2)5{'7^-~qi)~qi2-~vi25{qi2-(yi^h) (C12) 

XO (-gi2 • 1^12) Fi^h,l{V^;l {XI,X2) 

^ \ ^ I dxidx2~qi2 ■ 'V12S {qi2 - c^hh) © (-$12 • ^12) Fhh,i'A,i {xi,X2) 

lll2l[l'2 

xfhh {X1X2) [6 i/r - ~qi) -5{'r - "g's)] 

where use has been made of the the asymmetry of Fi^i^^ii^i'^-i (xi, X2) and of the symmetry of the 
distribution under an interchange of atoms. Finally, write 

- qi) - 5{r - q2) = / dx—5ir -xqi-il-x) q 2) (C13) 

Jo dx 

= -V • ~qi2 I dx 6(1^ - x~qi - {1 - x) ~q2) 
Jo 

so that the balance equation becomes 

j^^i , t) + V • ^ (r , t) (r , t) + v-Fi (r , t) = Si (r , t) (cu) 

with the flux written as = Ff('r^,t) + Fj{'r',t) where the kinetic contribution is 

t) = j d^i fi , ^1, t) Vi^i (^1) (C15) 

and the coUisional contribution is 

(~r, t) = - / dxidx2 ~qi2 Cq 12 ■ ^12) S {qu - ai^i^) 6 (-512 • ^"12) (C16) 

xfhh {X1X2) Fi^i2,i{i'2;i {xi,X2) / dxSiV - x~qi - (1 - x) ~q2) 

Jo 

and the source is 

5"; (7", = - ^ / dxidx2 Cq 12 ■ ^12) S {qu - cri^i^) © {-qi2 ■ ^12) (C17) 
xfhh {X1X2) Si^i^^i'^i'^.i {xi,X2) S{'r' - 'q'l) . 

2. Local number density 

Setting t/ji (^1) = 1 one has that 

Fhh,i'A;i {xi,X2) = ^ (^kIiIi {X12) {Sw^ - Sii'^) - 6i,i'^6i^i'^ {611^ - 611^)^ (C18) 
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From the normalization condition 

1 = E<'2 (C19) 



lb 

one has that 



J2 FwAi (^1' - ^ht (^12) {^w, - Sir, - Sih + h,) (C20) 



2 , , 

'1^2 '1^2 



2 , , 

^1^2 '1^2 



so that the balance equation becomes 



jni + V- (^nO + V • 7i = S-J") (C21) 

with the source 

^i"'^ — ~\'n J dxidx2 {'q'i2-~vi2)S{qi2-ai^i2)Q{-qi2-~vi2) (C22) 

ablil2 

xfhh (X1X2) 5(1^ - ~qi) K^i^ {X12) {Sal + hi - Sn^ - Su^) 
and the number current j i— j f + jj with 

7f = y d^i fi (r, ^1, t) V^. (C23) 

and 



J I 



E y (ia;i(ia;2 "^"12 ■ ^12) 5 {qi2 - ai^i^) 6 (-$12 ■ ^12) (C24) 

ablil2 

xfhh {X1X2) K^X (^12) {Sal - Sbi - Sii^ + Sii^) / dxS{'r^ - x~qi - {1 - x) V2) • 

Jo 

The balance equations for total number and mass density follow immediately. Summing over / 
gives 

jn (7^, + V • ^ , t) n , + V • J] 7f = (C25) 

since the sum of the coUisional contributions to the number current vanishes. Similarly, multiply- 
ing by mi and then summing gives the balance equation for local mass density 

j^p (r, t) + V • ^ (r, t)p{V,t) + V -Q = S^Pl (C26) 
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Here, the anomalous mass flux is 



2 



I1I2 o ' I1I2' 



Q = ^ y (ia;i(ia;2 "g^i2 ("5^12 • ^12) (gi2 - (J/iiz) © (-^12 • ^12) (C27) 

x/ziJa (a:;iX2) -fr,"^^ (^12) {ma -nib- +mi2) dx5{'r - x~qi - {1 - x) "^2) 

Jo 

and the mass source term is 

^i""^ = ^ XI y ^^I'^^s ('g'12 • ^12) S {qi2 - ai^i^) 6 {-qu ■ ^12) (C28) 

ablil2 

xfHi2 {X1X2) 6 (T^ - ti) (^-12) Smf^i^, 
which is only nonzero if the collisions do not conserve mass. 

3. Momentum density 

Taking (V'l) = rrii'v'i = p^i in Eq.(C14) and summing over / gives 

—pH + V ■ {pUli ) + V • + g ^) = 'S^^ (C29) 
with P — P ^ + P^ + P ^ where the kinetic contribution is 

"pK^J^^f ^ ^^^^ (C30) 

I 

To explicitly write the remaining flux and source terms, we need 

{b'^Pi - P 1) = ^ ( 7 - Sm^Vu) (C31) 

giving 

= -l^^ill {^12 - ^) - Ismlllt (C32) 



At this point, it is useful to separate 7 ^^^^ iiito two parts: its value in the local rest-frame of the 
colliding atoms, 7 i^i^, and the part coming from the Galilean transformation to the lab frame 
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so that we have that the coUisional part of the flux is 



~ ^ f dxidx2~q 12 Cq 12 -^12)3 {qi2-cri-,i2)Q{-qi2-lj' 12) (C33) 

lll2l[l'2 

xfhh ixiX2) kIiIi (xu) ~Ji\il [ dx5{'r - x~qi - {1 - x) ~q2) 



and the contribution from the instantaneous exchange of mass is 

-P^ = ^ / dxidx2 ~q 12 Cq 12 ■ ^ 12) S{qi2-(Ji^i2)Q{-qi2-^ 12) (C34) 

hl2l[l'2 

fhh (^1^2) kIiIi {X12) (^Vi2-^^ {mi>^ - mi^ - mi>^ - mi^) J dx 6 {V - x~qi - {1 - x) ~q2) 



X 



and the source can be written as 



5 (p) = TtS^"^ + S (C35) 
f(p) 1 



S = 2 / ^^1^^2 ( g 12 • 12) ^ {qi2 - CThh) (-$12 • V 12) 

x(y -H^ fi,i, {X1X2) Kilfl (xu) Smlfj {V - ~qi) 



By using the balance equation for total mass density 

^plt + V • {It pit) = p^li + pit ■ VI? - ^ V • Q + liS^"^ (C36) 



we can wnte 



^ + ^- V^ + p-M V • P + g • Vi?j . (C37) 



dt 

4. Energy 

Taking {l^i) = Imivf in Eq.(C14) and summing over I gives 



^E+V ■ {liE) + V • F = S^^^ (C38) 



(C39) 



where the kinetic part of the flux is 

= Y.\'^'j dirifi{V,i'i,t)Vivi 

= ^ / dlti fi (r, ^1, t) Vi (Vi + Ity (C40) 
I 

= ~q^ + li-y^ 
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with the kinetic contribution to the heat flux being defined as 



(C41) 



The source term comes from the even part of the collision kernel 



c'1'2 



+ 



1 



2m; 



2m, 



6E 



;' /' 1 

'l'2 _j_ _ 



2 '1'2 



1 



(C42) 



1 



2m/, 2mi2 



hh 



SB'S + - 



u 



so that 



1 ;/ ;/ 

bm}? 

2 '1*2 



f(p) 



— 



A hh 



2 



with the (rest-frame) source term 

i ^ \ ^ I dxidx2 {~q 12 ■ ^ 12) S {qi2 - (^hh) Q {-Q12 ■ ^ 12) 



l\l2l']l'2 



X 



;' ;' 1 ;' ;' / — > >\ • 



fhh {X1X2) kIiIi (X12) S{T - Vi) 



The flux comes from the odd part of the collision kernel 



77'''l'2 



2 V '1'^ 



1 2 1 

Pi 



- 1 



1 



m/i + mi, 



{pI-pD- 



mil - f^h 



1 



pI 



mil + \2mi^ 2mi^ 



-pI 



The first term gives 

- 1) 

_ 1 

+ m;^ 
1 

m^i + mi^ 



m/, + m;. 



Pi + 2^W2 



1 /' /' — > 

2 '1'2 



P2-i^l{ii2 



1 /' /' — > 
-5m 1^1/ 

2 '1'2 



- v v 
dm// 

'1'2 



mi'^ + m;^ m;, + mi^ 



{pI-p^) 



(C43) 



(C44) 



(C45) 



(C46) 
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while the second is 



_ / 1 1 A ^^^^ 



■mi{ - mi 



(/' /' 1 ;' ;' 
hh 2 

mil — mv ( YY 1 /'/' 
^ + -bm^^V 




so 



2f;7 = y 12 • ^7? + P^i2 ■Vi2 + — SEy + -Smylw.l (C48) 

\mi'^ + m;^ y mi^ \mi'^ + m;^ / mi^ 
which gives, after some algebra. 



2Fl^^ = ^12 • 7 + ' l^hiy (C49) 



This can also be written as 



M (Tf ^\ , mymi^-mymi^ 



+ 



-^^hil - 7; {-mi[ + my^ - mi^ + mi^) (Vi2 - 'u) 



m[ + my^ 

+ [my^ - mi^ - my^ + mi^) ((^12 ^) + ^ {my^ - mi^ - my^ + mi^](S\) 



■ '1*2 



+ • 7 h/2 



so that 



where the different pieces of the heat flux vector are the usual coUisional contribution 



(C52) 



"g"^ = J dxidx2~qi2{~qi2-^i2)5{qi2- (Ti^h)Qi{-qi2-^i2) (C53) 

xfhh {X1X2) kIiIi (xu) (^12 - ^) • dx6{'r - x~qi - {I - x) 2) , 
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a part arising from the instantaneous transfer of mass 



mi'^mi, - mi'^mi^ 



X J dxidx2 ~qi2 Cq 12 ■ ^12) S (gi2 - cr^^^J © (-$12 • ^^12) 

xfhh {X1X2) kIiI^ {X12) f^hhV^ 

X dx 5 ("r^ — x~qi — {1 — x) ~q2) , 
Jo 

and a part arising from the loss of energy 



(C54) 



dxidx2 ~qi2 {'q'12 ■ ^12) S {qi2 - (Thi^) O (-^12 ■ 'vi-if255) 



2 ^-^ mil + rriji 

xfhh {X1X2) kIiIi (xu) SeIiIi / dx6{'r' - x~qi - (1 - x) 2) ■ 

^0 

A little rearrangement allows us to write the energy balance equation as 

^E+V ■{liE) + V ■-q+'^ ■['u -"P^ + V ■ Q^^'q) ^^ + li-^^\^u^S^''\ (C56) 

Alternatively, noting the relation between the total energy and the kinetic temperature 

D , „ 1 



E = -nkeT + ^pu"" 



(C57) 



gives an equation for the evolution of the kinetic temperature 

: Vli +V ■~q 



dt 



e- (C58) 
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